library(rgdal)
library(leaflet)
library(sp)
library(dplyr)
#Unzip files
unzip("CO-wells.zip", exdir = ".")
list.files("CO-wells")
## [1] "CO-HUC.cpg" "CO-HUC.dbf" "CO-HUC.prj" "CO-HUC.qpj"
## [5] "CO-HUC.shp" "CO-HUC.shx" "CO-wells.csv"
#Read in data and convert to SpatialPointsDataFrame
CO.Wells <- read.csv(file ="CO-wells/CO-wells.csv", header = T, na.strings = "NA")
wells <- SpatialPointsDataFrame(coords = CO.Wells[,2:1], data = CO.Wells, proj4string = CRS("+proj=longlat +datum=WGS84"))
#Read in polygons and change change proj4string
huc1 <- readOGR("CO-wells", "CO-HUC")
## OGR data source with driver: ESRI Shapefile
## Source: "CO-wells", layer: "CO-HUC"
## with 93 features
## It has 2 fields
huc <- spTransform(huc1, CRS("+proj=longlat +datum=WGS84"))
#Records HUC ID for each point in wells
wells@data$HUC_ID <- over(wells, huc)$HUC_ID
#Groups by HUC ID and records total wells, active, and inactive
huc.wells <- wells@data %>% group_by(HUC_ID) %>% summarize(total = length(Facil_Stat), active = sum(Facil_Stat == "AC" | Facil_Stat == "PR"), inactive = total - active)
## Warning: package 'bindrcpp' was built under R version 3.4.2
huc.wells
## # A tibble: 91 x 4
## HUC_ID total active inactive
## <int> <int> <int> <int>
## 1 873 1 0 1
## 2 880 5 0 5
## 3 899 773 179 594
## 4 903 313 139 174
## 5 932 203 31 172
## 6 936 2 0 2
## 7 947 3095 1404 1691
## 8 951 4308 1945 2363
## 9 963 1231 571 660
## 10 964 4 0 4
## # ... with 81 more rows
#Joins huc.wells and the data of huc by HUC_ID
huc@data <- left_join(huc@data, huc.wells, by = "HUC_ID")
dim(huc@data)
## [1] 93 5
head(huc@data)
## HUC_ID HUC_NAME total active inactive
## 1 846 Upper North Platte NA NA NA
## 2 873 Upper Green-Flaming Gorge Reservoir 1 0 1
## 3 880 Upper Laramie 5 0 5
## 4 899 Little Snake 773 179 594
## 5 903 Vermilion 313 139 174
## 6 932 Lower Lodgepole 203 31 172
#Replace NA values with 0
huc@data[is.na(huc@data)] <- 0
huc@data
## HUC_ID HUC_NAME total active inactive
## 1 846 Upper North Platte 0 0 0
## 2 873 Upper Green-Flaming Gorge Reservoir 1 0 1
## 3 880 Upper Laramie 5 0 5
## 4 899 Little Snake 773 179 594
## 5 903 Vermilion 313 139 174
## 6 932 Lower Lodgepole 203 31 172
## 7 936 Upper Lodgepole 2 0 2
## 8 947 Crow 3095 1404 1691
## 9 951 Cache La Poudre 4308 1945 2363
## 10 963 Lone Tree-Owl 1231 571 660
## 11 964 Lower Green-Diamond 4 0 4
## 12 972 Lower South Platte 35 0 35
## 13 977 Sidney Draw 281 5 276
## 14 983 North Platte Headwaters 756 98 658
## 15 985 Upper Yampa 938 122 816
## 16 989 Lower Yampa 591 45 546
## 17 996 Middle South Platte-Sterling 6536 272 6264
## 18 1007 Pawnee 2970 646 2324
## 19 1011 Stinking Water 171 69 102
## 20 1022 Lower White 4601 1354 3247
## 21 1032 Upper White 522 134 388
## 22 1033 Big Thompson 2228 1142 1086
## 23 1034 Upper Republican 1 0 1
## 24 1036 Middle South Platte-Cherry Creek 23080 12248 10832
## 25 1038 Frenchman 391 176 215
## 26 1051 Colorado headwaters 77 0 77
## 27 1067 Piceance-Yellow 2239 795 1444
## 28 1068 North Fork Republican 6375 2682 3693
## 29 1086 St. Vrain 4944 2868 2076
## 30 1097 Kiowa 1188 197 991
## 31 1100 Beaver 3179 110 3069
## 32 1105 Bijou 1809 121 1688
## 33 1114 Blue 0 0 0
## 34 1124 Colorado headwaters-Plateau 12615 8328 4287
## 35 1128 Parachute-Roan 5575 3099 2476
## 36 1131 Eagle 2 0 2
## 37 1138 Clear 45 3 42
## 38 1141 Arikaree 2335 1086 1249
## 39 1149 South Fork Republican 800 313 487
## 40 1150 Westwater Canyon 10 0 10
## 41 1152 Upper South Platte 24 0 24
## 42 2207 Roaring Fork 33 10 23
## 43 1161 Little Beaver 40 0 40
## 44 1173 South Fork Beaver 71 1 70
## 45 1176 South Platte Headwaters 28 0 28
## 46 1182 Arkansas Headwaters 17 0 17
## 47 1187 North Fork Gunnison 187 30 157
## 48 1196 Lower Gunnison 222 8 214
## 49 1207 Big Sandy 1560 210 1350
## 50 1215 Lower Dolores 2 0 2
## 51 2208 East-Taylor 1 0 1
## 52 1225 Fountain 49 0 49
## 53 1226 Rush 355 27 328
## 54 1234 North Fork Smoky Hill 111 2 109
## 55 1236 Chico 56 0 56
## 56 1240 Upper Gunnison 12 0 12
## 57 1241 Uncompahange 48 0 48
## 58 1245 Horse 82 0 82
## 59 1250 Smoky Hill Headwaters 493 100 393
## 60 1258 Upper Arkansas 424 59 365
## 61 1264 Upper Dolores 205 10 195
## 62 1265 San Miguel 246 86 160
## 63 1266 Tomichi 1 0 1
## 64 1277 Upper Arkansas-Lake Meredith 51 0 51
## 65 1278 Ladder 440 100 340
## 66 1282 Upper Arkansas-John Martin 1053 86 967
## 67 1291 Whitewoman 239 19 220
## 68 1295 San Luis 22 0 22
## 69 1302 Saguache 3 0 3
## 70 1308 Montezuma 211 31 180
## 71 1319 Lower San Juan-Four Corners 136 4 132
## 72 1328 Middle Arkansas-Lake Mckinney 8 0 8
## 73 1331 Huerfano 371 37 334
## 74 1332 Animas 1765 1294 471
## 75 1336 Rio Grande headwaters 4 0 4
## 76 1341 Apishapa 426 241 185
## 77 1350 Alamosa-Trinchera 16 0 16
## 78 1352 Purgatoire 3531 2516 1015
## 79 1353 Mcelmo 366 95 271
## 80 1355 Upper San Juan 2452 1547 905
## 81 1361 Two Butte 64 1 63
## 82 1365 Piedra 90 38 52
## 83 1371 Mancos 463 10 453
## 84 1380 Bear 137 3 134
## 85 1381 Middle San Juan 761 338 423
## 86 1400 Conejos 4 0 4
## 87 1402 North Fork Cimarron 335 46 289
## 88 1415 Sand Arroyo 167 36 131
## 89 1423 Rio Chama 3 0 3
## 90 1426 Upper Cimarron 307 75 232
## 91 1432 Upper Rio Grande 2 0 2
## 92 1435 Cimarron headwaters 13 0 13
## 93 1444 Canadian headwaters 19 0 19
pal <- colorBin("YlOrRd", domain = huc$total, bins = 6, pretty = TRUE)
#Create custom labels
labels <- sprintf(
"<strong>%s</strong><br/>%s<br/>%s<br/>%s",
paste("HUC:",huc@data$HUC_NAME), paste("Total:",huc@data$total), paste("Active:",huc@data$active), paste("Inactive:",huc@data$inactive)
) %>% lapply(htmltools::HTML)
leaflet(width = "100%") %>%
# Overlay groups
addPolygons(data = huc, fillColor = ~pal(total), group = "Oil & gas wells", fillOpacity = .4, color = "white", opacity = 1, weight = 2, dashArray = "3", highlight = highlightOptions(
weight = 5,
color = "#666",
dashArray = "",
fillOpacity = 0.4,
bringToFront = TRUE),
label = labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto",
noHide = TRUE))%>%
# Base groups
addProviderTiles(providers$OpenStreetMap.Mapnik, group = "OSM (default)") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Imagery") %>% addProviderTiles(providers$Esri.WorldTopoMap,
group = "Topo") %>%
# Layers control
addLayersControl(overlayGroups = c("Oil & gas wells"), baseGroups = c("OSM (default)",
"Imagery", "Topo"), options = layersControlOptions(collapsed = FALSE)) %>%
#Add Legend
addLegend("bottomright", pal = pal, values = huc$total, title = "Total number of gas and oil<br>wells by USGS HUC")